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Abstract 



We critique a Pade analytic continuation method whereby a rational poly- 
nomial function is fit to a set of input points by means of a single matrix 
inversion. This procedure is accomplished to an extremely high accuracy us- 
ing a novel symbolic computation algorithm. As an example of this method 
in action, it is applied to the problem of determining the spectral function 
of a single-particle thermal Green's function known only at a finite number 
of Matsubara frequencies with two example self energies drawn from the T- 
matrix theory of the Hubbard model. We present a systematic analysis of the 
effects of error in the input points on the analytic continuation, and this leads 
us to propose a procedure to test quantitatively the reliability of the resulting 
continuation, thus eliminating the black magic label frequently attached to 
this procedure. 
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I. INTRODUCTION 



Analytic continuation arises in the many-body problem whenever real time dynamics 
are to be recovered from a response function calculated at non-zero temperatures in the 
Matsubara formalism. In that case, the function whose value is known only at a discrete set 
of points on the imaginary axis must be continued to the real axis. 

A general statement of the problem of interest in this paper is as follows: An analytic 
continuation of a function / defined on a subset A C C is a function that coincides with / 
on A and is analytic on a domain containing A. Usually, we are interested in the analytic 
continuation / with the largest such domain, for then / is the greatest analytic extension of 
/ to the complex plane. Since there exists no general prescription for finding / from /, there 
is no choice but to resort to approximate techniques. Currently, the state of the "art" is to 
interpolate between known points using fitting functions capable of reproducing the analytic 
structure of / in the complex plane§. A serious difficulty is that the analytic structure of / 
is not usually known a priori. 

A widely used technique is the Pade approximant method in which ratios of polynomials 
(or terminating continued fractions) are used as fitting functions. Several Pade schemes 
exist. The most common scheme, a recursive algorithm called Thiele's Reciprocal Difference 
Method!, was used by Vidberg and Serene! in the context of the Eliashberg equations. Yet, 
despite twenty years of widespread use, the Pade approximant method remains somewhat 
of an untested approach in that there is still no reliable, quantitative measure of the quality 
of a Pade result. The prevailing wisdom is that a Pade fit can be considered 'good' when 
the output function is stable with respect to the addition of more input points. The results 
of this work make it clear that such a criterion is insufficient. 

The various Pade schemes can be divided into two broad classes: (I) those which return 
the value of the continued function point by point in the complex plane (f(A),z) i— »■ f(z) 
and (II) those which yield the function itself f(A) i— > / by returning the polynomial (or 
continued fraction) coefficients. Thiele's method is class I, as are most numerical methods. 
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In this work we present a robust Pade scheme that is class II and propose a goodness-of-fit 
criterion based on the convergence of the polynomial coefficients to allowed values. One 
advantage of our approach is that we formulate the problem as a matrix equation, allowing 
us to make use of existing, highly efficient routines for matrix inversion. In contrast, a 
naively implemented recursion algorithm can lead to a severe propagation of error since 
repeated operations are performed on terms of very different orders of magnitude. 

Our paper is organized as follows. In the next section we review the formal aspects of 
thermal Green's functions to establish the definitions of the various functions that enter into 
this problem. In §111 we present the details of the Pade form that we will use, and state the 
algorithm that we use to solve for the Pade coefficients. This leads to the consideration of 
the accuracy required for such a calculation, thus necessitating the use of a high accuracy 
symbolic computation algorithm. This is presented in §IV, and then we display our numerical 
results for relevant test functions, including the statistical test that allows us to conclude 
whether or not a given analytic continuation is accurate. Finally, in §V we present our 
conclusions. 



First, we introduce the components of theories based on thermal Green's functions to 
establish the definitions of the various functions that enter into this problem. 

The one-particle propagator or Green's function can be formulated using real or imagi- 
nary time operators. In real time, the retarded Green's function 



describes how the system responds when a particle is added at time zero and removed at 
time t. Its imaginary time counterpart, the thermal Green's function 



II. GREEN'S FUNCTION FORMALISM 



G R (t) = -t({c(tt),c\0)})9(t) 



(1) 



G(t) 



(T[c(r)ct(0)]> , 



(2) 
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is not so clearly physically motivated. Its main advantages are its mathematical elegance 
and computational ease. Further, since it is defined in terms of the time ordering operator 
T, G(t) admits a diagrammatic expansion via Wick's theorem. Moreover, whereas the 
retarded Green's function G R (t) is aperiodic in t (it has a lone discontinuity at t — 0), the 
temperature Green's function is periodic in r with period 2/3. 

The two Green's functions have Fourier representations: the first a Fourier transform 



G R (t) = / due- lwt G K {uo) (3) 



2tt 

and the second, as a consequence of its periodicity, a Fourier series 



G{t) = \Y. e~ tm7TT/l3 Gm (4) 



odd m 



= I^ e — GK) (5) 

which, in Eq. ([5]), we have recast as a sum over the Matsubara frequencies {u n = (2n— 1)tt//3 : 
n G Z} of some new Fourier component G(u n ). 

The formal connection between the real and imaginary time formalisms is the following: 
There exists a unique function G : C i— > C with asymptotic form 

G» = (l/z)(l + C(l)/Imz) (6) 

which takes on the values of the Fourier components of the temperature Green's function at 
Matsubara points on the imaginary axis G(iuj n ) = G(u n ) and gives the Fourier transform 
of the retarded Green's function just above the real axis G(lj + i0 + ) = G R (u). That is, 
Eqs. (H) and (|5]) can be written as 



2tt 

and 



oo 



G R (t) = — I due- iwt G(to + iO + ) (7) 



DC 
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Clearly, all the information one can potentially extract from these functions is contained in 
G. 

The function G has several interesting properties. First, it is analytic everywhere in the 
complex plane with the exception of the real axis; this is a causality requirement. Second, 
the value of G in the upper and lower half planes is related by G(z*) = G(z)*, which is 
a statement of the time reversal symmetry between the retarded and advanced Green's 
functions. Its immediate consequence is that the imaginary part of G may be discontinuous 
across the real axis. It also implies that we need only know the function in either the upper 
or the lower half plane since the other is a conjugated reflection of the first. Third, G can 
be written as a Stieltjes/Hilbert transform 

G (2)= r^iM (9) 



where the spectral function, given by the magnitude of the discontinuity in G across the real 
axis, viz. 

A(u) = --lmG R (cu) 

7T 

= — -(G(u + iO + ) -G(u-iO + )), (10) 
is non-negative and normalized to unity 

/oo 
duA(u) = l. (11) 
-oo 

Typically, we are working in the Matsubara formalism and we calculate G{uj n ) from its 
self-energy (via G(u n )~ l = iu n — £ — S(c<j n )) which is in turn calculated from an approximate 
theory based on, e.g., a diagrammatic expansion of the propagator. From here the route to 
real time dynamics is somewhat circuitous: 

G(u n ) ^ G(z) Z G r (lo) G R (t) . (12) 

(1) The first step is to analytically continue from the Fourier components of the temperature 
Green's function to construct G. That this is possible, in principle, provided we know 



G(u> n ) = G(iu n ) for an infinite set of points including the point at infinity, was proved by 
Baym and Merminl. (2) Supposing that the analytic continuation to the upper half plane 
can be found, we merely evaluate it along the real axis (setting z = oj + i0 + ) to get G R (cu). 
(3) A Fourier transform then recovers the real time response function. 

In practice, however, we do not know the values of G(u n ) at an infinite number of points. 
Moreover, even if we did, the theorem of Baym and Mermin shows only the existence of a 
function G. There is no general method to perform the analytic continuation — hence the 
need for a procedure such as the Pade. 

III. PADE APPROXIMANTS 

The Pade method is based on the assumption that G can be written as a rational poly- 
nomial or terminating continued fraction. Since theories are most commonly specified by a 
choice of self-energy, the continued fraction form turns out to be the more useful, at least 
for investigating questions of a mathematical nature (e.g. analytic structure). In particular, 
we shall find it helpful to consider G (in the upper half plane) a continued fraction of Jacobi 
formil (J-frac). That is, 

G(z) = G (p+x) (z) = — ^ ^ ^- (13) 

z — e — z — e\— z — e r 

' (14) 



Z-£ - E(r)(*) 

where the X n and e n are complex constants. By comparison with Dyson's equation, Eq. ( |14| ) . 
we make the identification A 2 , = 1 and eo = £, where £ is just the free particle energy 
measured with respect to the chemical potential. Then, we find that £( r )(z) is itself a 
continued fraction 

\ 2 X 2 X 2 

%)^) = — — ■ (15) 

z — e\— z — e2— z — e r 

The justification for this continued fraction form is a theorem due to Wall and Wetzeli 
which assures us that a positive definite J-frac has a spectral representation with non- 
negative, integrable spectral weight and that it is analytic in the upper half complex plane 
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- all the properties we know G must have to be physically reasonable. By positive definite 
J-frac we mean a continued fraction in the form of Eq. ( |T3D satisfying Ime„ < and for 
which there exists a sequence of real numbers go, gx, ■ ■ ■ (0 < g n < 1) such that 

(ImA n ) 2 = (Ime n _i)(Ime n )(l - g n -i)g n ■ (16) 

There are two special cases worth mentioning. If the X n and e n are all real then the 
J-frac is positive definite and can be cast as a sum of simple poles@ 

Ejtf (17) 

n=l Z 

with real, distinct energies E n and positive residues R n > 0. The J-frac is also positive 
definite if the A„ are real and none of the e n sits in the upper half complex plane (Eq. (|1^) is 
satisfied by setting all g n = or 1), in which case the function is characterized by simple poles 
resting on or below the real axis. In the general case, all the continued fraction coefficients 
have the potential to be complex, with the exception of Aq = 1, Cq = £, and Xf. Since 
eo = £ has no imaginary part, Eq. ([16]) implies that the coefficient Xf must always be real 
and positive. 

It is clear that by observing the values of the A n , e n coefficients, one can learn a great deal 
about the analytic properties of Gr r+ iy For example, if some e n has a positive imaginary 
part (and no X m = for m < n) then GV r+ i) may have a pole in the upper half plane — such 
a function would be noncausal and have negative spectral weight. (In fact, it is through 
such considerations that we are led to propose a method for testing the accuracy of a given 
analytic continuation via a Pade.) 

Nonetheless, despite the usefulness of the continued fraction form, for computational 
purposes it is actually much easier to work with rational polynomials. Conveniently, every 
terminating continued fraction is equivalent to a rational polynomial. For instance, a J-frac 
with r stories, Eq. ([15]) say, can be written as the ratio 



of two polynomials P, Q defined recursively by the formulas 

P(„+i) = - e n )P {n) {z) - A^P (n _i)(z) (19a) 

Q(n+1) = (z ~ e n )Q(n)(z) - \ 2 n Q(n-x){z) (19b) 
(for n = 1, 2, 3, . . .) with base cases 

P (0) = 0, P (1) = A? (20a) 

Q ( o) = 1, Q(i) = 2 - ei . (20b) 
Writing out the leading order terms of P and Q 

P (r) (z) = Xjz^ 1 - \j(e 2 + e 3 + • • ■ + e r )z r ~ 2 + ■■■ (21a) 

Q( r )(z) = z r - (d + e 2 + ■ • • + e r )z r_1 + ■ • • (21b) 

makes it clear that the polynomial P is of order r — I'm z while the polynomial Q is of order 
r. (Accordingly, one refers to E( r ) in Eq. as a [r — 1/r] rational polynomial.) Moreover, 
it suggests that we write the self energy explicitly as a rational polynomial of the form 

£(„(,)= + . (22) 

gi + + • • ■ + qr-z r - 1 + z r 

It is straightforward to relate the old and new coefficients to one another via Eqs. (|T9|) and 
(0): e.g. \j = p r , ei = p r -i/p r - Qr, etc. 

The coefficients p n , g n can be determined by specifying the value of E( r ) at 2r points, 



viz., by solving the set of 2r linear equations 

{S( r) (io; n ) = E(a; n )} . (23) 

If we define the column vectors 



S 



p 
q 



pi 

Pr 
Ql 

Qr 



and er 



0" 2r .UCo> 2r 



(24) 



where a n = £(u; n ) are the known values of the self-energy at 2r Matsubara frequencies, and 
a matrix 



X 



1 iuj\ ■•• {iuJi) r 1 ~ a i '" ~ (J i(^ u} iY 1 
1 iuj 2 ■■■ («a; 2 ) r-1 -cr 2 ••• -<7 2 (2CL> 2 ) r_1 



1 zu; 2r . 



\r-l 



(25) 



equivalent to the system of equations given by Eq. (p3|), then the entire process of analytic 
continuation is reduced to a single matrix inversion 



P 

q 



X- l a 



(26) 



which provides the polynomial coefficients necessary to construct 



S (r)0) 





1 z z 2 


■■■ z r ~ 


1 


P 


1 


z z 2 ■ ■ 




q + z r 



(27) 



What we propose is that, having determined the p n , q n coefficients, we recover the A n , e n 
coefficients and then use the criteria provided by Wall and Wetzel's theorem to determine 
whether the matrix inversion produced a Gjy+i) with an acceptable analytic form. As a first 
step, we investigate what can be learned from A^, the first non-trivial J-frac coefficient. is 
equal to the sum of the residues of the poles in the self-energy and as such it gives the high 
frequency asymptotic behaviour of the self-energy via S( r )(^) ~ X 2 /z. A necessary condition 
for positive definiteness is that X 2 be real and positive. We shall see that the convergence of 
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Im to zero as a function of the number r of poles in the Pade fitting function can provide 
information on the quality of the fit and on the analytic structure of the true continuation 
G. 



IV. NUMERICAL RESULTS 



The procedure we have outlined in Sect. Ill is a specialization of the following general 
Pade procedure — such considerations are central to our statistical analysis of the quality 
of the fits provided by this method. 

Given a function / and a set A of 2r input points, we suppose that we can approximate 
the analytic continuation / by a [r — 1/r] rational polynomial /( r ), the coefficients of which 
are determined by solving the linear system of equations {/( r )(a) = /(a) : a G A}. This 
problem can be cast as a matrix inversion in which the kernel X has elements with ratios 
as large as 



Thus to reliably perform the inversion we need a numerical range ~ ( 2 , i.e. 2 log 10 £ decimal 
digits of numerical precision. This analysis is general in that no other Pade algorithm can 
have less stringent precision requirements. 

For the case of a self-energy £, known at the first 2r Matsubara frequencies above the 
real line on the imaginary axis, we have shown that the matrix X is given by Eq. fl25|) . Since 
E(ov) ~ 1/k-Vij t ne ratio of the largest to smallest terms in X is ( = {uJ2r) r — ((4r — l)7rT) r , 
the square of which gives an estimate of the amount of precision needed to invert X. Here, 
that corresponds to 



decimal digits. 

To achieve a sufficient level of precision for our numerical work, we implement the Pade 
algorithm using the symbolic computation package MAPLE. Under MAPLE, expression 




(28) 



2rlog 10 (4r - l)nT 



(29) 
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evaluation takes place in software and thereby transcends the limits imposed by hardware 
floating-point. All computations are performed in base ten to any desired level of precision 
(we specified Digits := 250 Moreover, MAPLE is an ideal environment for rapid 

prototyping since high level matrix data types and routines are available as primitives. 
We begin by considering a test function of known analytic structure. The self-energy0 

w„) = E "w)G (Q - k, v n , - u n ) (30) 

Q V rJ 



corresponds to the first 'rung' of the ladder diagrams in the T-matrixI3 approximation of 
the single-band Hubbard mo dell! (characterized by a near-neighbour hopping integral t and 
an on-site repulsion energy U). Here, G° is the free propagator 

G°(k,u n )= , 1 (31) 

and x° is the free pair susceptibility 

X°(Q, fn) = ^E G°(k, u n ,)G°(Q -k,u n - u n .) . (32) 

k a V 

The frequency sums in Eqs. (^) and (|32|) can be performed analytically, givingB 
and 

E( *"^ = ^g *4. + «o_ £ -fe-^ ' (34) 

Since the ^ are real, the analytic continuation of the self-energy is a meromorphic function 
with a finite number of simple poles, all situated along the real axis. Calculated in two 
dimensions on an 8 x 8 (M = 64) lattice, its k = component possesses r = 26 poles. 
(The number of poles is determined by counting the number of distinct elements in the set 
:V#,<3}.) 

For a particular set of parameters^ — we use an interaction strength \U\/t = 4, chemical 
potential fi/t = —2, and temperature T/t = 0.7 — the test function, Eq. (|30|), is calculated 
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in two different ways for the Matsubara frequencies {u>i,u>2, ■ ■ ■ ,u)2r}- First, it is calculated 
exactly, as prescribed by Eq. fl34|), but with a small, random error, viz., each value is mul- 
tiplied by 1 + e with — 1 < e < 1. Second, it is calculated by truncating the Matsubara 
sum at an arbitrary cutoff frequency v p ^> 1 (much larger than the relevant energy scale of 
the problem) and then systematically adding back the high frequency contributions up to a 
given order. That is, 

- ^ ^ X°(Q, Vn')G°(Q ~ k, Vn> ~ W n ) = ~ ^ ^ (Q , V n >)G (Q ~ k, V n > - UJ n ) 

v n> Wn'\< U P 

m—1 



+ E x? (Q)e (m) K + $5_ s ] + o(i/(v p ) m ) (35) 



where 



fe 

are the coefficients of a Laurent expansion of X°(Qi u n) and the 0(n functions (defined in 
the Appendix) are constructed using the symbolic manipulation capabilities of MAPLE. 

Now, we let the self-energy, evaluated at the first 2r Matsubara frequencies according to 
the two schemes described above, serve as the input to the Pade procedure. The resulting 
approximant E( r ) yields a propagator G( r+ i) (k, z) = (z — £g — S( r )(/c, z))^ 1 with spectral 
function Af r+ x\(k,u}) = — (l/7r)Im(5( r , +1 )(A;, u + iO + ). The spectral function derived from the 
Pade approximant is compared to that of the exact function using the logarithmic measure 



10- F = dx{A(k,x) - A (r+1) (k,x)) 2 

J — oo 

1 f°° / - -> - - \ 2 

= —J dx Im yG(k, x + irj) — G( r+ i)(k, x + irj)J . (37) 

In practice we choose rj to be a small, but noninfinitesimal positive real quantity (we use 
rj/t = 0.064), which has the effect of introducing a slight artificial broadening to the 5- 
function peaks of the spectral function. 

The results of this comparison (for the k = component of the spectral function) are 
presented in Figs. 1(a) and 1(b), where F is plotted as a function of r for different values of 
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the random error E = — log 10 e and the systematic error E = — log 10 l/(z/ p ) m = mlog 10 v p 
(and thus a larger E corresponds to a smaller error). In each graph, a vertical dashed line 
marks the exact number of poles (tq = 26) in the true self-energy. The most distinctive 
feature of both graphs is that, at high accuracy (large E), the F curves exhibit a large step 
at the point r = r . In the random error case, the E = 120 curve jumps by four decades, 
and this represents an improvement in the Pade fit of nearly 40 orders of magnitude. In the 
systematic error case, the result is even more dramatic: the E = 100 and E = 120 curves 
jump by roughly four and seven decades, respectively. 

At these large accuracies, the only factor inhibiting the success of the Pade approximants 
is the lack of a sufficient number of poles to reproduce the analytic structure of the true 
function. The large jump observed in the large E curves marks the point, r = r , at which 
the number of poles in the Pade approximant exactly matches the required number, and 
for this and larger r there is no difficulty in finding an excellent fit of the test function. 
In contrast, when the input points are known to relatively low accuracy, no such feature is 
observed, and instead the F curves pass smoothly through r . This makes clear that for 
self-energies calculated to 20, 40, or even 60 decimal digits of accuracy, the level of error in 
the input points is still the main obstacle to a successful Pade fit. 

The usual response to this situation is to increase the number of Pade points in an 
attempt to overcome the intrinsic error limitations (by making the system of equations more 
and more overcomplete) . However, whatever advantage this additional information brings 
to the Pade approximant is soon outweighed by the accompanying complications: When a 
rational polynomial of degree [r — 1/r] is used to fit a function with r < r poles, r — r zeros 
of the numerator must coincide with an equal number of zeros in the denominator in order to 
cancel the extraneous poles. As r — r grows, it is less and less likely that this cancellation 
will be complete. A slight misplacement of zeros leads to 'defects' in which the function 
moves between and oo in a small neighbourhood. Moreover, it cannot be predicted where 
these zero-zero pairs will appear^. For the purposes of calculating a spectral function, they 
are of little consequence provided that they lie deep in the complex plane. However, when 
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they are not so far removed from the real axis, they can distort the spectral function away 
from its proper shape. When they lie on or near the real axis, they can give rise to deep 
troughs of negative spectral weight and other spurious, non-physical features. 

The deterioration of the Pade fit, as described above, is evident in Figs. 1(a) and 1(b) 
in which many of the F curves reach maxima at points r best > r and then quickly begin to 
fall off for larger r. Interestingly, this behaviour is much more pronounced in the systematic 
error case where such maxima occur for each curve. In the random error case, the curves 
below some error threshold are essentially flat for all r. 

The primary lesson that one should draw from these results is that the addition of Pade 
points well beyond the required number is not a useful strategy for improving the Pade fit. 
Unless the exact analytic continuation is already known, there is no way to predict the value 
of r bcst- We believe that better results are achieved by fixing the number of Pade points 
at 2ro (giving rise to a [r — l/r ] rational polynomial) and working towards increasing the 
accuracy with which those input points are calculated. Even a small effort there can result 
in an improvement of several orders of magnitude in the fit. What to try when one does not 
know a priori what r is discussed later in this paper. 

Now consider Figs. 2(a) through 2(d) in which the spectral function of a Pade approx- 
imant with 26 poles (calculated by specifying the value of the self-energy at 52 Matsubara 
frequencies) is compared to the exact spectral function. In Fig. 2(a), the accuracy of the 
input points is given by E = 16 (random error), roughly the number of digits in a double 
precision Fortran variable. Despite the fact that the overall energy scale is correct, the 
details of the fit are quite poor. Here, the effect of insufficient accuracy is to produce a 
washed out version of the spectral function which completely lacks fine structure. Even at 
E = 30 (Fig. 2(b)), corresponding to the number of digits available in the largest Fortran 
data type, the Pade inversion is only just beginning to distinguish the main peaks of the 
spectral function. Figure 2(c) shows the result for E = 80 and Fig. 2(d) the result for 
E = 120. Notice that in Fig. 2(d), the fit is near perfect: even the smallest peaks have been 
reproduced faithfully. 
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In this example, with r = ro, the Pade approximant provides a remarkable fit to the 
true function whenever the accuracy of the input points is better than E ~ 110. The 
difficulty in translating our success in this specific case to the general problem is that, in 
real applications, one has no way to judge when sufficient accuracy has been achieved. Also, 
in most instances, the number of poles in the self-energy is unknown. 

In what follows, we hope to address these deficiencies. We begin by defining a logarithmic 
measure of the imaginary part of the J-frac coefficient A^: 

10" A =|ImA?|. (38) 

We argued in Sect. Ill that A^ ought to be real and positive. In a Pade calculation, however, it 
is real- valued only to within some small fraction which characterizes the numerical sensitivity 
of the matrix inversion. As we shall soon discover, the convergence of the imaginary part 
of A^ to zero (A — > oo) can be used (1) to determine when the threshold of accuracy for an 
exact fit has been reached and (2) to infer the value of r if it is unknown. 

In Figs. 3(a) and 3(b), we plot A as a function of r for the random and systematic error 
cases. Over each plot is superimposed a reference line given by Eq. (|2Tj|). What we observe 
is a set of A curves that initially follow the reference line but later fan out, spaced according 
to their E values. Our claim is that these curves provide the quantitative measure of success 
of the Pade approximant that has heretofore been lacking, the essential point being that 
the shape of the curves reveals the performance characteristics of the Pade inversion in the 
various r regimes. 

When < r < r , the accuracy of the Pade approximant is matrix inversion dominated 
and the behaviour of A is governed by A ~ 2rlog 10 (4r — l)nT. In this regime, the Pade 
approximant has too few poles to fit the true function and thus the matrix inversion must 
judiciously arrange the available poles (sometimes apportioning one pole to a region where 
there should be two or three) to give the best possible fit. In the opposite limit, r ^> ro, the 
accuracy of the Pade approximant is input-point error dominated. In this regime, there are 
more than enough poles to perform an exact fit, but the proper placement of those poles and 
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the determination of their residues is hampered by the finite accuracy to which the input 
points are known. We find this reflected in the A curves which, for large r, saturate at a 
value A ~ E (roughly). 

Most interesting, though, is the behaviour of A in the vicinity of r = ro where the A 
curves in Figs. 3(a) and 3(b) first cross the reference line. In those plots, we see that the A 
curves corresponding to small values of E closely follow the reference line (Eq. fl29|) ) until 
finite accuracy becomes a limiting factor. The curves then fall below the reference line and 
become more or less flat. As E is increased, the r coordinate at which a given A curve 
first deviates from the reference line moves to the right until (for some accuracy, Eq say) it 
coincides with ro- Here, there is a sudden change in behaviour: all A curves corresponding 
to accuracies E > E cross the reference line at r = r . Such a crossing signals that there 
are now both sufficient poles in the approximant and sufficient accuracy on the input points 
to fit E more or less exactly. We can verify this interpretation by appealing to Figs. 1 and 
2 which clearly show a large jump at r for precisely the same curves that demonstrate a 
crossing in Figs. 3(a) and 3(b). 

The results we have described are extremely general and do not dependent on the choice 
of test function. For example, we may replace Eq. fl5C|) with the full non-self-consistent 
T-matrix self-energy 

5 v> i + £/x°(Q,*v) 

Here, the frequency sums cannot be performed analytically^ and thus we do not have a 
closed form analytical expression for the self-energy. (Thus, this is more representative of 
the usual situation in which the Pade method might be applied.) In this case, we know only 
that its analytic continuation has a finite number of poles along the real axis (although we 
are able to predict analytically an upper bound for the number of poles). 

This self-energy can be calculated to high accuracy using the method of Eq. (|35|) with the 
Xq)(Q) replaced by the coefficients of the Laurent expansion of x°(Q; u n)J (1 + Ux°(Q, ^n))- 
That is, x m (Q) ^ X( 1} (Q), X°(2)(Q) ^ X( 2) (<5) - U X ° {1) (Q) 2 and so on according tc§ 
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(40) 

The Pade approximant method can then be applied to Eq. (^) calculated in this way. 
We find that the resulting plot of A vs. r is identical to that of Fig. 3(b) except that the 
crossing of the reference line at high accuracy now occurs at r = 156. This allows us to 
deduce that the function has tq = 156 poles, significantly more than the 26 poles of Eq. ( pip . 
(This is a consequence of the lifting of degeneracy in each Q component brought about by 
the renormalization 1/(1 + U(Q, v^))-) We also find that the approximant spectral function 
compares well with increasing accuracy of the input points to the numerically exact spectral 
functions as calculated (i) by a non-Pade method due to Marsiglio et alH (this non-Pade 
method is of limited application since it requires the self-energy to have a very specific form, 
but for those cases where it is applicable, it can outperform the Pade method), and (ii) by 
an exact partial fraction decomposition of the self energy^ that can be done to a very high 
accuracy (say 10~ 40 on all poles and residues). 

Finally, one interesting feature that could potentially be exploited is that for self-energy 
values calculated using the function expansion, the value of r which gives the maximum 
value of A roughly tracks rb es t (cf. Figs. 1(b) and 3(b)). 

V. CONCLUSIONS 

The Pade procedure is very sensitive to the numerical precision with which the matrix 
inversion is performed and to the intrinsic error on the input points. Sufficient precision 
is difficult to achieve in traditional computer languages (e.g. C, Fortran) and so, in many 
instances, it may be necessary to make use of a symbolic computation package capable 
of supporting very large precision data types. Likewise, sufficient accuracy is difficult to 
achieve without a sophisticated computational scheme (e.g. the function expansion) that 
goes beyond a simple truncation of the Matsubara frequency sums in the self-energy. The 
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required level of precision and accuracy depends on the temperature T, which controls the 
spacing of the Matsubara points, and on the pole count r . 

An insufficient level of accuracy leads to an approximant spectral function that lacks fine 
structural detail or, worse, one that exhibits spurious spikes or troughs of spectral weight. 
This poses a problem whenever we are interested in the presence of a specific feature in the 
spectral function (e.g. the onset of a normal state pseudogap). In that case, it is essential to 
have confidence in the quality of the Pade result. We must be convinced that the observed 
feature is robust and not merely a by-product of insufficient accuracy. 

We have argued that simply adding more Pade points cannot compensate for too large 
an error on the input points. While there is a small set of r values for which an increase in 
r improves the fit, there is no known criterion that indicates when to stop adding points. 
Without already knowing the exact result, one cannot distinguish between the regime where 
additional points improve the fit (r < rb cs t) and the regime where such points degrade it 
( r > r bcst)- Instead, we recommend the use of a Pade approximant function having the same 
number of poles as the function to be fit. The exact number of poles, when it is not known, 
can be determined from the crossing point in a A vs. r plot. The crossing also indicates that 
a sufficient level of numerical accuracy in the input points has been acheived. 

There are several caveats to the procedure we have outlined. (1) If the true Green's 
function has a branch cut along the real axis arising from trancendental functions then no 
A crossing will ever be observed, since a branch cut of that kind can only be represented by 
an infinity of poles (r = oo). (2) The self-energy of the Green's function we are trying to 
reproduce must have the correct asymptotic form and must be analytic in, say, the upper 
half of the complex plane; otherwise, the rational polynomial (or continued fraction) form 
of the approximant cannot reproduce its analytic structure. (3) The Pade method is often 
used to model a function that is smooth in some region of interest (well away from its 
poles) and such calculations are rarely performed with more than machine accuracy. Our 
numerical analysis of the Pade inversion, with its prediction of extremely high accuracy 
requirements, is not meant to invalidate these results. We have applied the Pade method 
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to the particularly difficult problem of reproducing the sharp peak structures characteristic 
of a spectral function whose Green's function has its poles along the real axis. In that case, 
the poles lie in the region of interest. The precision and accuracy requirements of the Pade 
inversion are greatly reduced if the poles of the Green's function lie deep in the complex 
plane. 

Finally, let us remember that the starting point for our new Pade approach was the 
realization that the convergence of the continued fraction coefficients to 'allowed' values 
can provide a criterion for judging the quality of a Pade approximant, even if the analytic 
structure of the function we are trying to fit is unknown. In Sect. IV, we demonstrated the 
utility of this idea using the Ai coefficient. However, we know that there is much addtional 
information that can be extracted from the remaining continued fraction coefficients. In 
future, perhaps our analysis can be extended to include ex, A 2 , e 2 , etc. 
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APPENDIX: 



In addition to the usual occupation functions 

1 ^ e iuJ " 0+ 
(3 iuj n - x ' eP x + 1 



f[ x ]= yl = (Ala) 



1 p WnO + 1 

= --V = (Alb) 

it is often convenient to define partial occupation functions. For example, the bose version 



of such a function looks like 



I p iv n Q + I , R N 



/? ^— ' ZZ/ n — X 27TZ V 27TZ 



where ip(z) = d\nT(z)/dz is the digamma functionill. This can be generalized to a m-order 
function (symmetric in its arguments) 

b[x 1 ,x 2 ,...,x m ] = -~'y]- — — (A3) 

p ^— ' %v n - xx iv n - x 2 %v n - x m 

which has the interesting property that it can be expressed (via partial fraction decomposi- 
tion) in terms of the (m — l)-order partial occupation function 



b[xi,X2,...,x m -2,x m -i]-b[xi,X2,...,x m -2,x m ] 



if %m—l 7^ ^T, 



b[x 1 ,x 2 ,...,x m ] = { _ Xm - X Xm (A4) 

-^b[x h x 2 ,...,x m - 2 ,y}\ y=Xri 



■gjjb[x u x 2 , . . -,x m -2,y]\ v=Xm otherwise. 

Equation ( |A2| ) serves to terminate the recursion. 

Furthermore, it is straightforward to show that for all / > 

\Vn\>V v 



-I y (-L L 



P ^ \{iv n ) l+1 W n - X (-il> n ) l+l -il> n - X 

fr[ o,o,„.,q , x] + (-l)'6[p,0 ,...,0, -ac] 

2+1 2+1 

1 <9< 



n^{^' y] + ( " 1) ^ [ " x ' y] }L=o (A5) 



20 



where, according to Eq. ( |A4|) , the two- argument function b[x,y] is related to b[x] by 

j[«,,] = MzM (A6) 
x-y 

provided x ^ y. The B functions provide a closed-form representation of the high-frequency 
asymptotics of a broad class of Matsubara sums. In particular, the sum 

- \ £ *°(& Vn>)G\Q - fc, zv - uj n ) (A7) 
can be separated into a finite sum over all low frequencies 

~ £ X Q {Q,Vn>)G Q {Q-k,v n ,-uj n ) (A8) 
and an infinite sum over the remaining frequencies 



+ £ X?o(Q)©(i+i)[^n + (A9) 



1=1 

where, in Eq. we have used the fact that the free susceptibility x° admits a Laurent 

expansion in the frequency variable 

= ^ + ^| + xkf + ... (A10) 

with Q-dependent coefficients 

4)(^) = ^ £(/&] + /fe-d - i)(e s + . (ah) 
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FIGURES 

FIG. 1. For various levels of (a) random and (b) systematic error, characterized roughly as 
W E (see the text for more details), the quality of the Pade fit as measured by F (see Eq. (|37|)) is 
plotted with respect to the number of poles in the Pade approximant (the solid lines are a guide 
to the eye). The k = self energy being studied is that of Eq. ( |30| ) where the parameters of the 
attractive Hubbard model (with t being the hopping energy), for an 8 x 8 square lattice, are a 
repulsive energy \U\/t = 4, a chemical potential n/t = —2, and a temperature T/t = 0.7. The 
vertical dashed line indicates the number of poles (tq = 26) in the true Green's function. In plot 
(a) , error bars (representing the standard deviation of the data points over a set of initial random 
seeds) are smaller than the symbols marking the data points and are not shown. In plot (b), the 
dotted line is the best linear fit through the maximum values of F. 

FIG. 2. The k = spectral function of the Pade approximant is compared to the exact spectral 
function for different levels of random error (10~ E ) on the initial input points. The parameters of 
the Hamiltonian and the self-energy being studied are the same as those of Fig. 1. 

FIG. 3. For various levels of (a) random and (b) systematic error (10~ E ), the parameter A is 
plotted with respect to the number of poles in the Pade approximant. The parameters are the same 
as those of Fig. 1. The vertical dashed line indicates the number of poles (tq = 26) in the true 
Green's function. The solid line originating in the lower left corner is given by 2r log 10(4r — 1)ttT. 
In plot (b), the dotted line is the best linear fit through the maximum values of A. The parameters 
of the Hamiltonian and the self-energy being studied are the same as those of Fig. 1. 
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